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Abstract 

We introduce a new approach to constructing networks with real- 
istic features. Our method, in spite of its conceptual simplicity (it has 
only two parameters) is capable of generating a wide variety of net- 
work types with prescribed statistical properties, e.g., with degree- or 
clustering coefficient distributions of various, very different forms. In 
turn, these graphs can be used to test hypotheses, or, as models of ac- 
tual data. The method is based on a mapping between suitably chosen 
singular measures defined on the unit square and sparse infinite net- 
works. Such a mapping has the great potential of allowing for graph 
theoretical results for a variety of network topologies. The main idea of 
our approach is to go to the infinite limit of the singular measure and 
the size of the corresponding graph simultaneously. A very unique 
feature of this construction is that the complexity of the generated 
network is increasing with the size. We present analytic expressions 
derived from the parameters of the - to be iterated- initial generat- 
ing measure for such major characteristics of graphs as their degree, 
clustering coefficient and assortativity coefficient distributions. The 
optimal parameters of the generating measure are determined from a 
simple simulated annealing process. Thus, the present work provides 
a tool for researchers from a variety of fields (such as biology, com- 
puter science, biology, or complex systems) enabling them to create a 
versatile model of their network data. 

1 Introduction 

As our methods of studying the features of our environment are becoming 
more and more sophisticated, we also learn to appreciate the complexity of 
the world surrounding us. The corresponding systems (including natural. 
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social and technological phenomena) are made of many units each having an 
important role from the suitable functioning of the whole. An increasingly 
popular way of grabbing the intricate structure behind such complex systems 
is a network or graph representation in which the nodes correspond to the 
units and the edges to the connections between the units of the original 
system [H [21 E]- It has turned out that networks corresponding to realistic 
systems can be highly non-trivial, characterized by a low average distance 
combined with a high average clustering coefficient anomalous degree 
distributions [5l [6] and an intricate modular structure [3 El [9]. A better 
understanding of these graphs is expected and, in many cases have been 
shown, to be efficient in designing and controlling complex systems ranging 
from power lines to disease networks [TU] . 

As increasingly complex graphs are considered a need for a better rep- 
resentation of the graphs themselves has arisen as well. Sophisticated visu- 
alization techniques emerged |Tl] and a series of new parameters have been 
introduced over the years [Il[2l|3]. Very recently one of us (L.L.) proved that 
in the infinite network size limit, a dense graph's adjacency matrix can be 
well represented by a continuous function W{x, y) on the unit square [121 113] ■ 
A similar approach was introduced by BoUobas et al. [HI [15] and used to 
obtain convergence and phase transition results for inhomogeneous random 
(including sparse) graphs. This two variables symmetric function (which can 
have a very simple form for a variety of interesting graphs, and was sup- 
posed to be either continuous or almost everywhere continuous) predicts the 
probability whether two nodes are connected or not. (The non-trivial rela- 
tions between the limiting objects of graph sequences and 2d functions are 
discussed in more details in the Supporting Information). In this paper we 
develop the above ideas further in order to obtain simple and analytically 
treatable models of random graphs with a level of complexity growing to- 
gether with their size. This is an important conceptual step acknowledging 
a rather natural expectation: the internal organization of larger networks is 
more complex than those of the smallest ones. (E.g., the social contacts in 
large universities are much more structured than in an elementary school, 
which is in part due to the underlying hierarchical organization of almost 
every large networks we know of.) 

In a sense, using a function to represent a network is very much like using 
a model to describe a network. Models in the context of networks have been 
playing a crucial role since they are ideal from the point of grabbing the 
simplest aspects of complex structures and thus, are extremely useful in un- 
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derstanding the underlying principles. Models are also very useful from the 
point of testing hypotheses about measured data. Indeed, many important 
and successful models have been proposed over the past 10 years to interpret 
the various aspects of real world networks. However, a considerable limita- 
tion of these models is that they typically explain a particular aspect of the 
network (clustering, a given degree distribution, etc), and for each new -to 
be explained- feature a new model had to be constructed. 

In the recent years, generating graphs with desired properties has at- 
tracted great interest. A few remarkable methods have been proposed, in- 
cluding the systematic approach for analyzing network topologies by Ma- 
hadevan et al., using the dK-sehes of probability distributions [16] . These 
distributions specify all degree correlations within d sized subgraphs of a 
given graph, with OK reproducing the average degree, IK the degree distri- 
bution, 2K the joint degree distribution, etc. Several methods for generating 
random graphs having a predefined finite dK-senes were also given in [T6] . 
(with typically d < 3). Most important of these techniques is based on 
rewiring of the links, as this turned out to be the only efficient tool in prac- 
tice. 

The concept of characterizing a network via the frequencies of given sub- 
graphs (forming a series with increasing size) is at the heart of the exponential 
random graph model as well [T71 [HI [19]. In this approach a possible sub- 
graph g (e.g., a pair of connected nodes, a wedge of a pair of links sharing 
a node, a triangle, etc.) is assigned a parameter rjg related to the frequency 
of the sub-graph, and the probability of a given network configuration is 
assumed to be proportional to exp(^ Vg^g)^ where Ug denotes the number of 
sub-graphs occurring in the network. The 1] parameters for a studied network 
are usually estimated using maximal likelihood techniques. 

The dK-series method and the exponential random graph model can be 
viewed as bottom-up approaches: in the first order approximation of the 
studied network we concentrate on the frequency of the most simple object 
(an edge), when this is reproduced correctly we move on to a slightly more 
complex sub-graph and so on. The series of sub-graphs from small/simple 
to large/complex are ordered into a sort of hierarchy. However, in a realistic 
scenario we stop in the above process at a relatively early stage, since on 
one hand most important properties of the networks are usually reproduced 
already, on the other hand including "higher order" sub-graphs becomes 
computationally very expensive. 

Hierarchy, self-similarity and fractality are very important concepts when 
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describing complex systems in nature and society, and turned out to be rel- 
evant in network theory as well [201 Ell [22] • Very recently two important 
network models have been introduced which are intrinsicly hierarchical, yet 
show general features. Avetisov et al. proposed in |23] the construction of 
random graphs having an adjacency matrix equivalent to a p-adic random- 
ized locally constant Parisi matrix, one of the key objects in the theory of 
spin glasses P^]. This symmetric matrix has a hierarchic structure, and its 
elements are Bernoulli distributed random variables (taking the value of 1 
with probability and the value with probability 1 — g^, where 7 counts 
the hierarchy levels). An interesting feature of this construction is that any 
sub-graph belonging to a specific hierarchy level 7 is equivalent to an Erdos- 
Renyi random graph pSj , nevertheless the overall degree distribution can be 
scale- free. 

The Kronecker-graph approach introduced by Leskovec et al. is cen- 
tered around hierarchic adjacency matrices as well, however in this case the 
self similar structure is achieved by Kronecker multiplication as follows [26] . 
Starting from a small adjacency matrix A^, (where A\- = 1 if nodes i and 
j are linked, otherwise A\- = 0), at every iteration we replace each current 
matrix element by multiplied by the matrix element itself, hence enlarg- 
ing the matrix by a factor given by the size of A^. In the stochastic version 
of this model the elements of A^ are replaced by real numbers between 
and 1, and at the final stage of the multiplication process we draw a link 
for each pair of nodes with a probability given by corresponding element 
in the obtained stochastic adjacency matrix. According to the results, the 
Kronecker-graphs obtained in this approach can mimic several properties of 
real networks (heavy tails in the degree distribution, and in the eigenvalue 
spectra, small diameter, densification power law) simultaneously. Further- 
more, in [27] Leskovec and Faloutsos presented a scalable method for fitting 
real networks with Kronecker graphs. 

We note that link probability matrices similar to the previous examples 
can be also used for community detection as pointed out by Nepusz et al. 
in [2H1 121] • In their approach (inspired by Szemeredi's regularity lemma 
[30j ) the diagonal elements of the matrix give the the link density inside the 
corresponding communities, whereas the off-diagonal elements correspond to 
the link probabilities between the groups. 

In summary, a plausible classification of the emerging graph generating 
procedures/approaches involves the following types. Generating graphs as 
i) stochastic growth processes (e.g., [IJ), ii) as a process of connecting or 
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rewiring nodes according to prescribed probabilities ([U ED 1321 E21 El]), iii) 
accepting varying configurations with a prescribed probability [T71 [TSl 
[19], iv) by deterministically or stochastically obtaining its adjacency matrix 
from simpler initial matrix, [23], |26],[27], v) from a function W{x,y) on the 
unit square providing a value for the probabilities of node pair connections 
ini[I31[I3. 

Rewiring and the related construction techniques do not provide a clue 
how a complex network emerges from a simple rule. On the other hand, gen- 
erating a graph from a fixed function/measure does not result in networks 
with increasing complexity. Our approach can be considered as a combina- 
tion of iv) and v) (thus, combining their advantages), assuming that in the 
infinitely large network limit the right representation is a singular measure 
(nowhere continuous function). 

Thus, here we introduce a new method to constructing random graphs 
inheriting features from real networks. The main idea of our approach is to 
replace W{x, y) by a fractal (singular) measure (also called multifractal), and 
go to the limit of infinitely fine resolution of the measure and the infinitely 
large size of the generated graph simultaneously. Consequently, the complex- 
ity of the obtained network is increasing with the size. Another advantage 
of this approach is that the statistical features characterizing the network 
topology, e.g., the degree distribution, clustering coefficient, degree correla- 
tions, etc, can be simply calculated anal5d;ically. For generating networks 
with a given prescribed statistical feature (e.g., a given degree distribution), 
the optimal parameters of the generating measure defining the multifractal 
can be determined from a simple simulated annealing process. 

2 Model 

The network generation has three main stages in our approach: we start 
by defining a generating measure on the unit square, next we transform the 
generating measure through a couple of iterations into a link probability mea- 
sure, and finally, we draw links between the nodes using the link probability 
measure. The generating measure is defined as follows. We identically divide 
both X and y axis of the unit square to m (not necessarily equal) intervals, 
splitting it to w? rectangles, and assign a probability to each rectangle 
(i,j G [1,"^] denote the row and column indices). The probabilities must 
be normalized, YliVij = 1 symmetric Pij = Pji. Next, the link probabil- 
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ity measure is obtained by recursively multiplying each rectangle with the 
generating measure k times (which is equivalent to taking the k-th tensorial 
product of the generating measure). This results in m^^ rectangles, each as- 
sociated with a linking probability Pij{k) equivalent to a product of k factors 
from the original generating pij. In our convention k = 1 stands for the 
generating measure, thus, a link probability measure at = 1 is equivalent 
to the generating measure itself. Finally, we distribute N points indepen- 
dently, uniformly at random on the [0, 1] interval, and link each pair with a 
probability given by the Pijik) at the given coordinates. The above process 
of network generation is illustrated in Fig{T], whereas in Figj2l we show a 
small network obtained with this method. 

We note that our construction could be made even more general by re- 
placing the "standard" multifractal with the fc-th tensorial product of a 
symmetric 2d function < W{x,y) < 1 defined on the unit square. Al- 
though the resulting Wk{xi, Xfc, yi, yt) = W{xi, yi) ■ ■ ■ W{xk, yk) func- 
tion is [0, l]^'^ — )■ [0, 1] instead of [0, 1]^ — )■ [0, 1], with the help of a measure 
preserving bijection between [0, 1] and [0, l]'^ it could be used to generate 
random graphs in the same manner as with our multifractal. 

The diversity of the linking probabilities Pij{k) (and correspondingly, the 
complexity of the generated graph) is increasing with the number of itera- 
tions, just like in case of a standard multifractal. In order to keep the gen- 
erated networks sparse, we must ensure that the average degree, {d) of the 
nodes does not change between subsequent iterations. This can be achieved 
by an appropriate choice of the number of nodes as a function of k, using the 
following relation: 

I- h 

m m 

i=l j=l 

where aij{k) denotes the area of the box i,j at iteration k. In the special 
case of equal sized boxes aij{k) = m"^*^, and due to the normalization of the 
linking probabilities the above expression simplifies to (d) = Nmr"^^ . Thus, 
to keep the average degree constant when increasing the number of iterations 
for a given generating measure, the number of nodes have to be increased 
exponentially with k. 
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Figure 1: Schematic illustration of the multifractal graph generator, a) The 
construction of the link probability measure. We start from a symmetric 
generating measure on the unit square defined by a set of probabilities pij = 
Pji associated to m x m rectangles (shown on the left). In the example shown 
here m = 2, the length of the intervals defining the rectangles is given by /i 
and I2 respectively, and the magnitude of the of the probabilities is indicated 
by both the height and the color of the corresponding boxes. The generating 
measure is iterated by recursively multiplying each box with the generating 
measure itself as shown in the middle and on the right, yielding m'' x m*^ 
boxes at iteration k. The variance of the height of the boxes (corresponding 
to the probabilities associated to the rectangles) becomes larger at each step, 
producing a surface which is getting rougher and rougher, meanwhile the 
symmetry and the self similar nature of the multifractal is preserved, b) 
Drawing linking probabilities from the obtained measure. We assign random 
coordinates in the unit interval to the nodes in the graph, and link each 
node pair J, J with a probability given by the probability measure at the 
corresponding coordinates. 
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Figure 2: A small network generated with the multifractal network generator, 
a) The generating measure (on the left) and the link probability measure (on 
the right). The generating measure consists of 3x3 rectangles for which 
the magnitude of the associated probabihties is indicated by the color in 
a similar fashion to FigH] The number of iterations, k, is set to k = 3, 
thus the final link probability measure consists of 27x27 boxes, as shown in 
the right panel. (Note that k = 1 corresponds to the generating measure 
in our convention), b) A network with 500 nodes generated from the link 
probability measure. The colors of the nodes were chosen as follows. Each 
row in the final linking probability measure was assigned a different color, 
and the nodes were colored according to their position in the link probability 
measure. (Thus, nodes falling into the same row have the same color). 
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3 Statistical Methods 



One of the main advantages of our model is that the statistical properties 
characterizing the network topology can be calculated analytically. An im- 
portant observation concerning our model is that nodes having coordinates 
falling into the same row (column) of the link probability measure are sta- 
tistically identical. This means that e.g., the expected degree or clustering 
coefficient of the nodes in a given row is the same. Consequently, the distri- 
butions related to the topology are composed of sub-distributions associated 
with the individual rows. 

Let us concentrate on the degree distribution first, which can be expressed 

as 



where p\ \d) denotes the sub-distribution of the nodes in row i, and Uik) 
corresponds to the width of the row (giving the ratio of nodes in row i 
compared to the number of total nodes). These pf'' can be calculated using 
the generating function formalism as shown in the Appendix, resulting in 



i. Even though the degree distribution of nodes in a given row follows a 
Poisson-distribution according to ([3]), the overall degree distribution of the 
generated graph can show non-trivial features, as will be demonstrated later. 

Similarly to the degree distribution, the clustering coefficient and the 
average nearest neighbors degree can be calculated analytically as well in 
a rather simple way (as given in the Supporting Information). According 
to Figl3]D-d, the analytical results for the quantities above are in very good 
agreement with the empirical distributions, (obtained by generating a num- 
ber of sample graphs for the chosen parameters). The use of analytic formulas 
instead of empirical distributions can significantly speed up the optimization 
of the generating measure with respect to some prescribed target property. 




(2) 



i=l 
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Figure 3: Comparison between the analytical and empirical result for a 
randomly chosen generating measure, a) The generating measure (left) and 
the link probability measure (right). The number of rows in the generating 
measure was set to m = 3 with equal box lengths U = 1/3, the corresponding 
initial linking probabilities pij were chosen randomly. The number of itera- 
tions was set to k = 4, resulting in the linking probability measure shown on 
the right. This link probability measure was used to generated 100 samples 
of random networks with N = 5000 nodes each, b) The degree distribution 
obtained by averaging over the samples (symbols), plotted together with the 
analytical result obtained from eq.Q (continuous line), showing very good 
agreement. The error-bars (showing the standard error of the mean) are 
smaller than the symbols. The average clustering coefficient (C) of nodes 
falling into the same row of the final link probability measure and the near- 
est neighbors average degree as a function of the node degree are plotted in 
panels (c) and (d) respectively. 
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4 Results 



Depending on the choice of the generating measure and the box boundaries, 
our method is capable of producing graphs with diverse properties. However, 
to generate a random graph with prescribed features in our approach we need 
to optimize the generating measure with respect to the given requirements. 
Let us suppose that the number of nodes in the graph to be generated is 
given. In this case we have two parameters: the number of boxes in the 
generating measure (given by m^), and the number of iterations, k. The 
actual Pij and box boundaries are "self adjusting", as we shall describe in 
the following. 

Let us denote the property to which we are optimizing the generating 
measure by J-'. A conceptually simple example is when our goal is to obtain 
a network with a given degree distribution, in this case is equivalent to 
p{d). In principle, T depends on p^j, 1^, k, N (and in an implicit way on m, 
through the box sizes and linking probabilities). However, as m, k and are 
kept constant, we discard them from the notation and write the "value" of the 
property corresponding to a given choice of p^j and k as J^{pij, h). (Note that 
in most cases J^{pij,li) is actually a high dimensional object, e.g, a degree 
distribution, and not a real number). The target value of the property to 
which we would like the system to converge is denoted by J-'*. 

In order to be able to make the studied property of the generated network 
converge to the goal J-'*, we have to define a way to judge the quality of 
the actual J^{pij, h). In other words, we have to define a sort of distance or 
similarity between T{pij, k) and T*. This distance/similarity measure can be 
used as an energy function during a so called simulated annealing procedure, 
and we shall denote it by E[J-'{pij,li),T*]. The actual form of this function 
depends on the actual choice of the property, e.g., in case of optimizing the 
degree distribution a plausible choice is the sum of the relative differences 
between the degree distributions: 

where d runs over the degrees, p'^^\d) is the value of the actual degree dis- 
tribution at degree d and p*{d) is the value of the target degree distribution 
at the same degree. 

In the simulating annealing we also define a temperature, T, which is de- 
creased slowly during the process. The process itself consist of many Monte- 
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Carlo steps, and in one step we try to change one of the hnking probabihties 
or one of the box boundaries by a small amount, following the Metropolis al- 
gorithm [35] . If the energy E2 after the change is smaller than the energy Ei 
before, the change is accepted. In the opposite case, the change is accepted 
by a probability given by P = exp[— (£'2 — Ei)/T]. 

The above procedure can be generalized in principle to optimizing with 
respect to multiple properties simultaneously as well. However, for simplicity 
here we consider the optimization of the different properties separately. In 
FigJl^ we show the results for optimizing the generating measure with respect 
to various target degree distributions. Although the three chosen targets 
are rather different, (a scale-free distribution, a log-normal one, and a bi- 
modal distribution), our method succeeded in finding a setting of Pij and k 
producing a degree distribution sufficiently close to the target. Similarly, in 
FigSb the results from optimizing with respect to three clustering coefficient 
distributions are displayed showing again a reasonable agreement between 
the targets and the results. 

5 Discussion 

Our approach raises a number of fundamental graph theoretical and prac- 
tical questions. Should we expect that large real graphs converge to some 
limiting network in a strict sense of the convergence? Or, alternatively, their 
structure cannot be mapped onto a fixed function, and only an ever changing 
(with the size of the network) measure (in the infinite network size limit be- 
coming singular) can be used to reflect the underlying structural complexity? 
This would be in contrast with the consequences of the renowned Szemeredi 
Lemma [30] valid for arbitrary dense graphs. 

Although it can be shown analytically (see SI) that in the infinitely large 
network size limit our construction converges to a relatively simple graph, the 
convergence to this structure is extremely slow. According to our numerical 
studies, there is a very extensive region between the small and infinite regimes 
in which a well defined, increasingly complex structure emerges as our method 
is applied. Details about aspects of the slowness of convergence involving 
an extremely slow growth of the relative number of isolated nodes and the 
appearance of oscillations are given in the SI. 

In summary, our results demonstrate that it is possible to use simple 
models to construct large graphs with arbitrary distributions of their essential 
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Figure 4: Optimizing the generating measure with respect to different target 
properties. During the optimization process the number of nodes, N, the 
number of rows in the generating measure, m, and the number of iterations, 
k, are kept constant, only the probabihties Pij and the length of the intervals 
li defining the generating measure are adjusted. The typical value of the 
constant parameters in our experiments were = 10000 — 20000, m = 3 — 4 
and k = 3 — 5. a) Optimizing with respect to different degree distributions. 
The target distributions are shown with circles, whereas the corresponding 
results at the end of the optimization procedure are marked by squares. 
The black symbols come from an experiment where the target was a power- 
law degree distribution (the inset shows this on log-scale), the gray symbols 
correspond to a setting with a log-normal target, whereas the white symbols 
show the results of an experiment with a bi-modal target distribution, b) 
Optimizing with respect to different clustering coefficient distributions. 

characteristics, such as degree distribution, clustering coefficient distribution 
or assort at ivity. In turn, these graphs can be used to test hypotheses, or, 
as models of actual data. The combination of the tensorial product of a 
simple generating measure and simulated annealing technique leads to small 
(in practice 3x3 to 5x5) matrices representing the most relevant statistical 
features of observed networks. A very unique feature of this construction is 
that the complexity of the generated network is increasing with the size. In 
addition, the multifractal measure we propose is likely to result in networks 
displaying aspects of self-similarity in the spirit of the related findings by 
Song et al. [2]. 
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Appendix 

Al Thermodynamic limit of graph sequences 

The infinitely large system size limit is very important in statistical physics, 
as the properties of the studied phenomena are manifested in the most pure 
way in this limit. Similarly, in complex network theory the infinitely large 
limit graph of a converging graph sequence can be considered as a "platonic" 
network exhibiting the fundamental common properties of the graphs in the 
sequence in the most pure way. But under what conditions can we say that a 
given sequence of graphs is converging to something non-trivial and yet suffi- 
ciently universal to be conceptually meaningful? A simple intuitive condition 
is that the statistical features commonly used to characterize a network (e.g., 
degree distribution, clustering coefficient, etc.) should converge. The actual 
definition is based on homomorphisms (adjacency preserving maps) as fol- 
lows. For two simple graphs F and G, let hom(F, G) denote the number of 
homomorphisms from V{F) (the nodes in F) to V{G) (the nodes in G). The 
homomorphism density t{F, G) is defined as the probability that a random 
map from V{F) to V{G) is a homomorphism, given by 



A sequence (Gn) of graphs is convergent., if the sequence t(F, G„) has a limit 
for every simple graph F. Losely speaking, this condition can be interpreted 
as the convergence of the probability of finding any given finite sub-graph in 
the sequence of networks. 

We note that when G has bounded average degree and F is connected, 
it is more natural to normalize the homomorphism density by dividing with 
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instead of Therefore, we let inj(F, G) denote the number 

of injective homomorphisms from graph F to graph G and define 



s{F,G) 



inj(F,G) 

\V{G)\ 



(6) 



(Thus, s{F,G) is the average number of labeled copies of F, such that a 
specified node of F goes on a specified node of G.) 

Convergent graph sequences are related to 2d functions in a non-trivial 
way fl2\ [T3] . First of all, we can construct a convergent graph sequence 
using a symmetric measurable function, < W{x,y) < 1, defined on the 
unit square as follows. For a given network size N, we distribute N points 
independently, uniformly at random on the [0, 1] interval. These points cor- 
respond to the nodes in the network, and each pair of nodes is linked with 
the probability given by W{x,y) at the coordinate of the according points. 
In the iV — )■ cxD limit the obtained graph sequence is converging. What is 
even more surprising, it can be proven that we can represent any convergent 
graph sequence by a 2d function, since for any convergent graph sequence 
one can find a W{x,y) providing the same limiting sub-graph densities. 

The average degree of nodes in a random graph generated from a given 
W{x,y) using the construction above can be given simply as 



Thus, in the N ^ oo limit the obtained network becomes dense. In con- 
trast, real networks are usually sparse in the sense that their average degree 
is not expected to grow with increasing size. A solution to this problem 
was proposed by Bollobas et al., by redefining the linking probabilities as 
W{x, y)/N, resulting in a network with an average degree independent of N. 
They showed that depending on the choice of W{x, y), a wide range of sparse 
networks can be generated. 

Our approach is different from this method in that instead of using a con- 
struction into which we build in the level of complexity from the beginning, 
we generate complexity by using tensorial products of increasing power as 
N CO. This is a qualitatively new picture, corresponding to reality to a 
higher degree (larger graphs are more complex/inhomogeneous/structured in 
nature than smaller graphs). In addition, we achieve this using a relatively 
simple construction. 
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Al.l Limiting cases of the multifractal model 

A shortcoming of our model is that it can lead to a network in which the 
majority of nodes are isolated in the N ^ oo limit. However, as we shall 
see, this effect is negligable for graphs in the size range of real networks. 

In general, if Wi{x, y), W2{x, y), . . . , Wk{x, y) is a sequence of symmetric 
measurable functions on the unit square (with < Wk{x,y) < 1 for any k), 
let us define Wk{x) as the average linking probability for a node at position 
X given by 

Wk{x) = / Wk{x,y)dy. (8) 

^0 

Similarly, let denote the average link probability for the whole network, 
which can be expressed as 



•^fc = / Wk(x)dx. (9) 
Jo 

Let us choose the number of nodes, Nk associated to Wk{x, y) in such a way 
that the average degree of nodes converges to a constant (non zero) {d) for 
k — >■ oo, thus 

NkUk ^ {d) . (10) 

(This means that the number of links is around {d) /2.) The degree dis- 
tribution of a node at position x can be given by a binoimal distribution 
as 

p{d, x) = {^f^ Wk{xY [1 - Wkix)]""^-'-" . (11) 

In the thermodynamic limit this can be approximated by a Poisson-distribution 
written as 

p(d,x)^ ^^^^;,^^^]' e-^^-^(-). (12) 

The degree distribution of the whole network is obtained by integrating 
p{d,x), resulting in 

p{d) = ^ J^' [NkWk{x)fe-''^^-^^'>dx. (13) 

In particular, the probability that a randomly chosen node will be isolated 
(having degree zero) is 



p{d = 0) = / e-^'=^'=(^)(ix. (14) 
Jo 



19 



From and (fTUj) it follows that the average value of Wk{x) is around (d) /Nj.. 
In case Wk{x) is actually independent of x, then Wk{x) = {d) /N^, and 

p{d = 0) e-^'^l (15) 

However, if Wk{x) is such that its typical value is much smaller than its 
average, then typically e~^'='"'=^^^ ~ 1 resulting in 

p(d = 0)~l, (16) 

which means that the majority of nodes becomes isolated. The condition for 
avoiding this degeneracy can be formulated as 

f e-^*""=(")cia; < c, (17) 
Jo 

where c < 1 is a constant. 

In case of the multifractal graph generator, (or a more general "tensoring" 
construction), the above condition is not fulfilled, unless Wk{x) is independent 
oix. As mentionned in the main text, by using a measure preserving bijection 
between [0, 1] and [0, l]'^, our model can be formulated in a more general form 
using the tensorial product = W^'' = W ■ ■ ■ W defined as 

W{xi, ...,Xk,yi,...,yk) = W{xi, yi) ■ ■ ■ W{xk, yk)- (18) 
The marginals ^ in this representation are given by 



Wk{xi,. . . ,Xk) = / W{xi,yi)---W{xk,yk)dyi...dyk 

J [0,1]'= 

= w{xi) ■ ■ -wixk), (19) 
where w{x) = W{x,y)dy. Similarly, is transformed into 

uJk= Wkixi, . . . ,Xk)dxi . . .dxk = u'', (20) 

J[0,1]'= 

where u = w{x)dx. Thus, according to (ITU]) , we should choose Nk — 

{d) /CJ^ 

Unfortunately, these functions do not satisfy condition f|T71) unless w{x) is 
constant. Indeed, if (xi, . . . ,Xk,yi, ■ ■ ■ , yk) is a random point in [0, l]'^, then 

\nwk{xi, . . . , Xk) = \nw{xi) + ■ ■ ■ + \nw{xk) ^ k / \nw{x)dx (21) 
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almost surely by the Law of Large Numbers. Let a = exp(J lnw{x) dx), then 
a < u by the Jensen inequality (expect if w is constant), and the value of 
Wk is almost always close to a'^, while its average is u''. Since (a/w)'^ — )■ if 
k ^ oo, this shows that if f lT7|) holds, then w{x) is constant. 

On the other hand if w{x) = u for any x, then the expected degree of the 
nodes becomes independent from their position and the degree distribution 
converges to a Poisson distribution, just like in case of an Erdos-Renyi graph. 
In this case it is also easy to calculate the number of copies of any connected 
graph F with / nodes in a graph Gk obtained from Wk- There are Nk{Nk — 
1) ■ ■ ■ {Nk — I + 1) ^ Nl ways to map V{F) into [uk] injectively, and for 
each map, the probability that it is a homomorphism is t{F, Wk) = t{F, W)^. 
Hence ^ 

E(inj(F,G.))~iV^t(F,Ty)^ = (ci)'(^t^) . (22) 

Since we have a sparse graph, we want to normalize this by A^^; so the 
normalized number of copies of F is 

E(.(F. G.)) = '221^ . ^--HiF, Wf = (d)'-' (^)' , (23) 

For example, the normalized number of triangles is 

l.(K,,a,^^,,f{^S!^^)\ (24, 

It is easy to see that if w{x) = u is constant, then 

t{F,W)<J~\ (25) 

where equality holds if and only if F is a tree or W is an equivalence relation 
such that there is a partition U • ■ ■ U Sm = [0, 1] and 

W{x,y) = l^^ ^^y^S. for some ^, ^^6) 
lO, otherwise. 

From fl23|) and fl25|) we gain 



{dy ^ , if F is a tree or W is an equivalence relation, 
0, otherwise. 

(27) 



EisiF,Gk))^ 
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Using high concentration inequahties one can prove that this convergence 
happens almost surely (not just in expectation). We see from (P7|) that the 
sequence Gk is convergent with probability 1 in the Benjamini-Schramm 
sense. 

There are a number of possibilities which one could try to cure the degen- 
eracy of the thermodynamic limit of our model shown here, however these 
are out of the scope of the present study. We could modify the tensoring con- 
struction by adding to W®'^ a constant Ck tending to reasonably slowly. 
Another possibility is to modify W®'' to (VT®*^)'^*, where ak — >■ 0. 

A1.2 Numerical studies 

Next, let us investigate the magnitude of the above effect for graphs in the 
size range of real networks. For this purpose we generated networks from 
randomly chosen generating measures (with m = 4, equal sized boxes) it- 
erated from A; = ltoA;=ll. The number of nodes at k = 1 was set to 
1000 and to 5000 respectively, and for /c > 1 it was adjusted using eq. () in 
the main text. (Thus the average degree of the graphs remained the same 
during the iterations.) In FiglSl we show the results obtained by averaging 
over 1000 samples for both settings by plotting p{d = 0) as a function of N. 
Inspite of the increasing tendency of the curves, at the last iteration with 
network sizes above 10^, the fraction of isolated nodes is still very low. Thus, 
the effect of isolated nodes becoming dominant is negligible on the scale of 
real world applications. 

In spite of the analytical results for the convergence in the thermodynamic 
limit, the degree distribution often shows an oscillatory behavior in the size 
range of real networks. This is shown in Figj6] for an m = 3 generating 
measure iterated from /c = 1, = 100 to A; = 11, = 4.24 ■ 10^. As k is 
getting larger, the more oscillations can be observed in p{d) towards the large 
degrees. 

A2 Calculating statistical distributions 

A serious advantage of our model is that the statistical properties character- 
izing the network topology can be calculated analytically. In the Appendix 
of the main text we give a derivation for the degree distribution, based on 
the generating function formalism. The definition and the most important 
properties of the generating functions can be summarized as follows. 
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Figure 5: The ratio of isolated nodes in function of the network size for 
graphs obtained from iterating random generating measures from A; = 1 to 
k — 11, averaged over 1000 samples. The number of nodes at iteration k = 1 
was chosen to be N{k = 1) = 1000 (circles) and N{k = 1) = 5000 (triangles). 
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Figure 6: The degree distribution obtained from Eq.(2) in the main text for 
an m = 3 generating measure iterated from k — 1,N — 100 to A; = 11, = 
4.24 • 10^. 



A2.1 The generating functions 

If a random variable ^ can take non-negative integer values according to some 
probability distribution — n) = p{n), then the corresponding generating 
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function is given by 



(28) 



n=0 



The generating-function of a properly normalized distribution is absolute 
convergent for all |a;| < 1 and hence has no singularities in this region. For 
X = 1 it is simply 



(29) 



n=0 



The original probability distribution and its moments can be obtained from 
the generating-function as 



, , 1 d'^GJx) 



oc 



a;=0 



n=0 



(30) 
(31) 



x=l 



And finally, if = + C2 + ••• +6) where Ci,C2, are independent random 
variables (with non-negative integer values), then the generating function 
corresponding to Virj — n) = (T{n) is given by 



G,{x) = (x") = ( n^^O = n (^^'> = G,,{x)G,,{x) ■ . ■ G,,{x). (32) 



\i=l 



i=l 



A2.2 The degree distribution 

The degree distribution of the nodes falling in row i of the link probability 
measure, pf\d), can be calculated as follows. In our construction we draw 
links for a node in row i pointing to nodes in row j altogether nj{k) times 
with a probability Pij{k), where nj{k) is the number of nodes in row j, given 

by 

rijik) = Nljik). (33) 

Since the link draws are independent, the distribution of the number of links 
from a node in row i to nodes in row j is binomial. This can be approximated 
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by a Poisson-distribution when Uj is sufficiently large as 



where (dij) denotes the average number of links from a node in row i to nodes 
in row j given by 

{dij) = nj{k)pij{k). (35) 

The degree of a node in row i is given by the sum over the links towards the 
other rows as 

di{k) = J2d^J{k). (36) 

j 

Therefore, the generating function of pf\d) is the product of the generating 
functions of the plj\d) distributions: 

G['\x) = llG[fix), (37) 



where G[j\x) is defined as 



oo oo 
d d=0 



{dijjk))'^ ^_(d^^(k))^d 



(38) 



(A summary of the most important properties of the generating functions 
is given in the Supporting Information). By substituting (!38|) into (137|) we 
arrive to 

j 

where we used that due to the independence of the links, the expected degree 
of a node in row i can be expressed as 

(d^) = E (d^,) ■ (40) 
j 

25 



An alternative form for (rfj) can be obtained by substituting fl35p into the 
equation above, yielding 



{d,) = NJ2p^Jik)^k). 



(41) 



j 



The degree distribution of the nodes falling into row i can be obtained by 
transforming back the generating function in fl39p . resulting in 



A2.3 The clustering coefRcient 

Similarly to the degree distribution, the clustering coefficient of nodes falling 
into the same row of the link probability measure is expected to be the same. 
The clustering coefficient of a node in row i can be obtained by calculating 
the number of triangles containing the node, divided by the number of link 
pairs originating from the node. Since the triangles are equivalent to link 
pairs originating from the node having their other end connected by a third 
link, the expected clustering coefficient of a node in row i can be given as 



The first term in both the numerator and the denominator corresponds to 
the link pairs for which the other end of the links point to the same row j, 
whereas the second terms give the contribution from link pairs connecting 
our node in row i to distinct rows j and q. 

A 2. 4 The average nearest neighbors degree 

Finally, we mention that the degree correlations can be calculated bom pij{k) 
(and li{k)) as well. Here we derive the expression for the average nearest 
neighbors degree, d^N, as a function of the node-degree. This is one of the 
most simplest quantity characterizing the degree correlations: an increasing 
curve corresponds to an assortative network, whereas a decreasing one signals 
disassortative behavior. The average degree of the neighbors of a node from 




{d.{k))' 



e 



{d,{k)) 



(42) 



(43) 
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row i can be given as 

The average degree of the neighbors of a node with degree d can be given 
as a sum over the possible (iyvTVi' muhiphed by the conditional probability 
p'^^\i\d) that the node is from row given that its degree is d: 

dPM = Y.P"'^('\d)dNN,. (45) 

i=l 

These conditional probabilities can be obtained as follows. The number of 
nodes from row i with degree d is ni{k) p[''\d) , whereas the total number of 
nodes with degree d is np'^^\d). The probability that a node is from row i 
given that its degree is d is the ratio of these two: 

By substituting and dS]) into (|15]) we get 

.w,., 1 E,"L\p.(fc)^.(fc)(^.(fc)) 

dNNW = A / ^i(^)Pi ('^) — ^ k • (47) 



27 



